%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Reviews of Accelerator Science and Technology
%% Trim Size: 11in x 8.5in
%% Text Area: 9.25in (include runningheads) x 6.6in
%% Main Text: 10/13pt
%% Date:      16-04-2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\documentclass[wsdraft]{ws-rast} % to draw visible frames for the text
%\documentclass[twocolumn]{ws-rast}
\documentclass[]{report}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{url}
\usepackage{hyperref}

\usepackage[displaymath]{lineno}

\input{newcommands}

\begin{document}

\markboth{J.-L. Vay, R. Lehe}{Simulations for plasma and laser acceleration.}

\title{WarpX}

%\author{Jean-Luc Vay, R\'emi Lehe}
%\address{Lawrence Berkeley National Laboratory, Berkeley, California, USA\\
%\email{jlvay@lbl.gov}}

%\address{Lawrence Berkeley National Laboratory, Berkeley, California, USA\\
%\email{rlehe@lbl.gov}}

\maketitle

\linenumbers

\begin{abstract}
Computer simulations have had a profound impact on the design and understanding of past and present plasma acceleration experiments, and will be a key component for turning plasma accelerators from a promising technology into a mainstream scientific tool. In this chapter, we present an overview of the numerical techniques used with the most popular approaches to model plasma-based accelerators: electromagnetic Particle-In-Cell, Quasi-Static, Ponderomotive Guiding Center.
The material that is presented is intended to serve as an introduction to the basics of those approaches, and to advances (some of them very recent) that have pushed the state-of-the-art, such as optimal Lorentz boosted frame, advanced laser envelope solvers and the elimination of numerical Cherenkov instability. The Particle-In-Cell method, which has broader interest and is more standardized, is presented in more depth. Additional topics that are cross-cutting such as azimuthal Fourier decomposition or filtering are also discussed, as well as potential challenges and remedies in the initialization of simulations and output of data.  Examples of simulations using the techniques that are presented have been left out of this chapter for conciseness, and because simulation results are best understood when presented together - and contrasted with - theoretical and/or experimental results, as in other chapters of this volume.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The electromagnetic Particle-In-Cell method}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}
%\begin{centering}
\includegraphics[scale=0.6]{PIC.png}
%\par\end{centering}
\caption{\label{fig:PIC} The Particle-In-Cell (PIC) method follows the evolution of a collection of charged macro-particles (positively charged in blue on the left plot, negatively charged in red) that evolve self-consistently with their electromagnetic (or electrostatic) fields. The core PIC algorithm involves four operations at each time step: 1) evolve the velocity and position of the particles using the Newton-Lorentz equations, 2) deposit the charge and/or current densities through interpolation from the particles distributions onto the grid, 3) evolve Maxwell's wave equations (for electromagnetic) or solve Poisson's equation (for electrostatic) on the grid, 4) interpolate the fields from the grid onto the particles for the next particle push. Additional ``add-ons'' operations are inserted between these core operations to account for additional physics (e.g. absorption/emission of particles, addition of external forces to account for accelerator focusing or accelerating component) or numerical effects (e.g. smoothing/filtering of the charge/current densities and/or fields on the grid).}
\end{figure}

In the electromagnetic Particle-In-Cell method \cite{Birdsalllangdon},
the electromagnetic fields are solved on a grid, usually using Maxwell's
equations

\begin{subequations}
\begin{eqnarray}
\frac{\mathbf{\partial B}}{\partial t} & = & -\nabla\times\mathbf{E}\label{Eq:Faraday-1}\\
\frac{\mathbf{\partial E}}{\partial t} & = & \nabla\times\mathbf{B}-\mathbf{J}\label{Eq:Ampere-1}\\
\nabla\cdot\mathbf{E} & = & \rho\label{Eq:Gauss-1}\\
\nabla\cdot\mathbf{B} & = & 0\label{Eq:divb-1}
\end{eqnarray}
\end{subequations}
given here in natural units ($\epsilon_0=\mu_0=c=1$), where $t$ is time, $\mathbf{E}$ and
$\mathbf{B}$ are the electric and magnetic field components, and
$\rho$ and $\mathbf{J}$ are the charge and current densities. The
charged particles are advanced in time using the Newton-Lorentz equations
of motion
\begin{subequations}
\begin{align}
\frac{d\mathbf{x}}{dt}= & \mathbf{v},\label{Eq:Lorentz_x-1}\\
\frac{d\left(\gamma\mathbf{v}\right)}{dt}= & \frac{q}{m}\left(\mathbf{E}+\mathbf{v}\times\mathbf{B}\right),\label{Eq:Lorentz_v-1}
\end{align}
\end{subequations}
where $m$, $q$, $\mathbf{x}$, $\mathbf{v}$ and $\gamma=1/\sqrt{1-v^{2}}$
 are respectively the mass, charge, position, velocity and relativistic
factor of the particle given in natural units ($c=1$). The charge and current densities are interpolated
on the grid from the particles' positions and velocities, while the
electric and magnetic field components are interpolated from the grid
to the particles' positions for the velocity update.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Particle push}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A centered finite-difference discretization of the Newton-Lorentz
equations of motion is given by
\begin{subequations}
\begin{align}
\frac{\mathbf{x}^{i+1}-\mathbf{x}^{i}}{\Delta t}= & \mathbf{v}^{i+1/2},\label{Eq:leapfrog_x}\\
\frac{\gamma^{i+1/2}\mathbf{v}^{i+1/2}-\gamma^{i-1/2}\mathbf{v}^{i-1/2}}{\Delta t}= & \frac{q}{m}\left(\mathbf{E}^{i}+\mathbf{\bar{v}}^{i}\times\mathbf{B}^{i}\right).\label{Eq:leapfrog_v}
\end{align}
\end{subequations}
In order to close the system, $\bar{\mathbf{v}}^{i}$ must be
expressed as a function of the other quantities. The two implementations that have become the most popular are presented below.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Boris relativistic velocity rotation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{Particle_pushers/Boris_pusher.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Vay Lorentz-invariant formulation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{Particle_pushers/Vay_pusher.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Field solve}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Various methods are available for solving Maxwell's equations on a
grid, based on finite-differences, finite-volume, finite-element,
spectral, or other discretization techniques that apply most commonly
on single structured or unstructured meshes and less commonly on multiblock
multiresolution grid structures. In this chapter, we summarize the widespread
second order finite-difference time-domain (FDTD) algorithm, its extension
to non-standard finite-differences as well as the pseudo-spectral
analytical time-domain (PSATD) and pseudo-spectral time-domain (PSTD)
algorithms. Extension to multiresolution (or mesh refinement) PIC
is described in, e.g. \cite{VayCSD12,Vaycpc04}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Finite-Difference Time-Domain (FDTD)}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{Maxwell_solvers/Maxwell_FDTD_solver.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Non-Standard Finite-Difference Time-Domain (NSFDTD)}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{Maxwell_solvers/Maxwell_NSFDTD_solver.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Pseudo Spectral Analytical Time Domain (PSATD)}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{Maxwell_solvers/Maxwell_PSATD_solver.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Current deposition}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{Deposition/Current_deposition.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Field gather}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input{Gather/Field_gather.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Filtering}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

It is common practice to apply digital filtering to the charge or
current density in Particle-In-Cell simulations as a complement or
an alternative to using higher order splines \cite{Birdsalllangdon}.
A commonly used filter in PIC simulations is the three points filter
$\phi_{j}^{f}=\alpha\phi_{j}+\left(1-\alpha\right)\left(\phi_{j-1}+\phi_{j+1}\right)/2$
where $\phi^{f}$ is the filtered quantity. This filter is called
a bilinear filter when $\alpha=0.5$. Assuming $\phi=e^{jkx}$ and
$\phi^{f}=g\left(\alpha,k\right)e^{jkx}$, the filter gain $g$ is
given as a function of the filtering coefficient $\alpha$ and
the wavenumber $k$ by $g\left(\alpha,k\right)=\alpha+\left(1-\alpha\right)\cos\left(k\Delta x\right)\approx1-\left(1-\alpha\right)\frac{\left(k\Delta x\right)^{2}}{2}+O\left(k^{4}\right)$.
The total attenuation $G$ for $n$ successive applications of filters
of coefficients $\alpha_{1}$...$\alpha_{n}$ is given by $G=\prod_{i=1}^{n}g\left(\alpha_{i},k\right)\approx1-\left(n-\sum_{i=1}^{n}\alpha_{i}\right)\frac{\left(k\Delta x\right)^{2}}{2}+O\left(k^{4}\right)$.
A sharper cutoff in $k$ space is provided by using $\alpha_{n}=n-\sum_{i=1}^{n-1}\alpha_{i}$,
so that $G\approx1+O\left(k^{4}\right)$. Such step is called a ``compensation''
step \cite{Birdsalllangdon}. For the bilinear filter ($\alpha=1/2$),
the compensation factor is $\alpha_{c}=2-1/2=3/2$. For a succession
of $n$ applications of the bilinear factor, it is $\alpha_{c}=n/2+1$.

It is sometimes necessary to filter on a relatively wide band of wavelength,
necessitating the application of a large number of passes of the bilinear
filter or on the use of filters acting on many points. The former
can become very intensive computationally while the latter is problematic
for parallel computations using domain decomposition, as the footprint
of the filter may eventually surpass the size of subdomains. A workaround
is to use a combination of filters of limited footprint. A solution
based on the combination of three point filters with various strides
was proposed in \cite{Vayjcp2011} and operates as follows.

The bilinear filter provides complete suppression of the signal at
the grid Nyquist wavelength (twice the grid cell size). Suppression
of the signal at integer multiples of the Nyquist wavelength can be
obtained by using a stride $s$ in the filter $\phi_{j}^{f}=\alpha\phi_{j}+\left(1-\alpha\right)\left(\phi_{j-s}+\phi_{j+s}\right)/2$
for which the gain is given by $g\left(\alpha,k\right)=\alpha+\left(1-\alpha\right)\cos\left(sk\Delta x\right)\approx1-\left(1-\alpha\right)\frac{\left(sk\Delta x\right)^{2}}{2}+O\left(k^{4}\right)$.
For a given stride, the gain is given by the gain of the bilinear
filter shifted in k space, with the pole $g=0$ shifted from the wavelength
$\lambda=2/\Delta x$ to $\lambda=2s/\Delta x$, with additional poles,
as given by $sk\Delta x=\arccos\left(\frac{\alpha}{\alpha-1}\right)\pmod{2\pi}$.
The resulting filter is pass band between the poles, but since the
poles are spread at different integer values in k space, a wide band
low pass filter can be constructed by combining filters using different
strides. As shown in \cite{Vayjcp2011}, the successive application
of 4-passes + compensation of filters with strides 1, 2 and 4 has
a nearly equivalent fall-off in gain as 80 passes + compensation of
a bilinear filter. Yet, the strided filter solution needs only 15
passes of a three-point filter, compared to 81 passes for an equivalent
n-pass bilinear filter, yielding a gain of 5.4 in number of operations
in favor of the combination of filters with stride. The width of the
filter with stride 4 extends only on 9 points, compared to 81 points
for a single pass equivalent filter, hence giving a gain of 9 in compactness
for the stride filters combination in comparison to the single-pass
filter with large stencil, resulting in more favorable scaling with the number
of computational cores for parallel calculations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Porting onto new architectures, parallelization, vectorization, mesh refinement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Just a reminder that you may have to run bibtex
%% All of it up to \end{document} can be removed
%% if you don't like the warning.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\IfFileExists{\jobname.bbl}{} {\typeout{} \typeout{{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}}
\typeout{{*}{*} Please run \textquotedbl{}bibtex \jobname\textquotedbl{}
to optain} \typeout{{*}{*} the bibliography and then re-run LaTeX}
\typeout{{*}{*} twice to fix the references!} \typeout{{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}{*}}
\typeout{} }



\end{document}
